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Abstract 

A  computer  model  is  developed,  capable  of  predicting  the  flow  distribution  along  the  height  of  a  fuel-cell  stack.  The  model  distinguishes  a 
number  of  hydraulic  resistances  that  are  linked  in  series  and  parallel,  thus  forming  a  network  to  simulate  the  gas  flow  and  pressure  distribution 
in  a  stack.  The  hydraulic  resistances  are  found  from  analytical  solutions  or  from  tabulated  data  inferred  from  measurements.  The  electrochem¬ 
istry  of  the  stack  is  simulated  by  adding  or  subtracting  gas  in  the  active  cell  area.  The  model  is  tested  on  an  assumed  separator-plate  geometry 
and  stacking  height. 
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i.  Introduction 

Fuel  cells  are  an  exciting  way  of  producing  electricity  by 
the  electrochemical  conversion  of  gases,  such  as  natural  gas 
or  hydrogen,  into  gaseous  products,  waste  heat  and  electricity. 
Because  the  reaction  takes  place  by  the  electrochemical  con¬ 
version  of  the  reactants,  and  not  via  a  thermal  reaction,  the 
efficiency  is  higher  than  that  produced  by  conventional 
means.  The  efficiency  can  be  further  increased  by  using  the 
waste  heat  in  co-generation  applications,  see  Refs.  1 1,2]. 

Solid  oxide  fuel  cells  (SOFCs)  are  fabricated  out  of  all 
solid-state  ceramic  components  and  have  many  advantages 
over  other  types  of  fuel-cell  systems,  including  higher  tem¬ 
peratures  for  co-geneiation  and  combined  heat  and  power 
(CHP)  applications,  all  solid-state  design  and  the  possibility 
of  using  a  wide  range  ofoxidizable  gases  (for  internal  reform¬ 
ing  processes)  [3,41. 

A  variety  of  configurations  has  been  considered  for  the 
&CFC  system,  including  the  planar,  monolithic  and  tubular 
designs,  although  there  are  many  other  variations  on  these 
three  basic  designs  [5]. 

Single  cells  are  usually  stacked  into  multi-cellular  units,  or 
stacks,  using  an  interconnect  or  bi  polar  plate.  The  form  of 
the  interconnect  plate  is  dependent  upon  the  design  of  the 
overall  fuel-cell  system.  In  the  planar  design,  for  example, 
the  interconnect  plate  can  be  used  to  support  the  single  cell, 
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to  act  as  a  bi-poiar  plate,  and  to  act  as  a  gas  manifold.  The 
primary  challenge  in  the  design  of  an  SOFC  stack  is  to  ensure 
that  the  interconnector  performs  the  tasks  required  of  it;  the 
major  design  criteria  is  the  distribution  of  the  gases  across 
the  electrode  faces  [6]. 

It  is  well  known  that  non-uniform  gas  distribution  will  have 
a  detrimental  effect  on  the  efficiency  and  performance  of  the 
SOFC  stack,  and  thus  optimization  of  the  gas-distribution 
profile  will  consequently  allow  for  a  higher  power  output 
[71. 

We  have  shown  [  8  ]  that  the  gas-flow  distribution  in  a  fuel¬ 
cell  stack  can  be  presented  in  an  analytical  form,  assuming 
that  the  stack  is  viewed  as  a  network  of  hydraulic  resistances 
to  (he  flow  of  gas.  Using  an  internally-manifolded,  planar 
stack  design  with  an  assumed  geometry,  the  ratio  between 
the  average  flow  and  the  flow  in  the  upper  cell  of  the  stack 
was  determined,  which  could  be  used  to  ascertain  separator 
plate  design. 

This  paper  looks  a(  a  numerical  approach  to  the  problem 
of  non-uniform  gas  distribution  in  an  internally-manifolded 
fuel-cell  stack.  The  feed  gases  and  waste  gases  are  assumed 
to  flow  in  opposite  directions,  which  based  Gn  current  tech¬ 
nical  insights  is  the  most  viable  option.  The  numerical  model 
is  two  dimensional,  as  it  describes  a  cross  section  running 
through  the  cells,  from  top  to  bottom,  including  the  feed  and 
waste  gas  channels.  The  novel  feature  about  this  model  is  that 
the  flow  distribution  is  considered;  most  of  the  modelling 
until  now  has  only  focused  on  the  modelling  of  planar  or 


216 


R.J  Boersma.  N.M.  Sammes  /Journal  of  Power  Sources  63  (1996)  215-219 


tubular  single  cells.  The  model  is  basically  a  numerical  ver¬ 
sion  of  the  analytical  model  described  in  Ref.  18],  but  it 
includes  a  larger  number  of  resistances  and  their  flow  depend¬ 
encies.  The  numerical  approach  has  the  advantage  that  a  large 
number  of  physical  phenomena  can  be  included.  Neverthe¬ 
less,  the  results  are  difficult  to  verify.  Therefore,  the  analytical 
approach,  as  presented  in  our  previous  paper  [8],  forms  a 
good  basis  to  verify  the  numerical  results. 


Z  Description  of  the  model 

To  calculate  the  flow  and  pressure  distribution,  the  stack 
is  modelled  as  a  series  and  parallel  connection  of  hydraulic 
resistances.  Each  separator  plate  is  modelled  as  consisting  of 
12  distinguishable  resistances,  linked  together  to  model  the 
stack,  as  shown  in  Fig.  1 . 

The  hydraulic  resistances  being  distinguished  are: 

R1  Resistance  to  flow  through  the  manifold  channel. 
This  is  modelled  as  a  flow  through  a  pipe  with  a 
rough  surface. 

R2  Resistance  attributable  to  the  splitting  of  gas  in  the 
direction  of  the  mainstream. 

R3  Resistance  attributable  to  the  splitting  of  gas  in  the 
direction  of  the  cell. 

R4  Resistance  to  the  flow  through  the  gas-distribution 
area,  which  is  between  the  actual  cell  and  the  mani¬ 
fold  channel. 

R5  Resistance  caused  by  the  contraction  of  the  flow 
when  entering  the  channels  adjacent  to  the  active 
area. 

R6,  R7  Resistance  caused  by  the  flow  of  the  gas  through  the 
channels  adjacent  to  the  active  area.  The  electro- 


(b) 

Fig.  1.  The  separator  plate  is  modelled  as  a  series  and  parallel  linking  of 
hydraulic  resistances.  The  model  distinguishes  12  resistances.  R1  through 
R12,  assigned  to  one  separator  plate.  Plan  (a)  shows  the  separator  plate 
from  the  top.  with  the  resistances  R3  through  RIO  projected  on  it.  Plan  <b) 
shows  the  cross  section,  through  the  manifold  channels,  of  the  plate,  thus 
identifying  the  resistances  R1 .  R2.  R1 1  and  R12. 


chemical  reaction  causes  the  gas  composition  to 
change,  and  thus  its  flow,  density  and  viscosity 
change.  This  effect  is  modelled  by  assigning  the  inlet 
properties  to  the  first  half  of  the  channel  and  the  outlet 
properties  to  the  second  half  of  the  channel.  The 
outlet  properties  are  calculated  on  the  basis  of  the 
mass  balance  and  the  thermal  balance,  which  follow 
from  the  current  through  the  cell  and  the  heat  pro¬ 
duced  in  the  cell. 

R8  Resistance  caused  by  the  expansion  of  the  flow  when 
leaving  the  channels. 

R9  Resistance  to  the  flow  through  the  gas-collection 
area,  which  is  between  the  actual  cell  and  the  outlet 
manifold  channel. 

RIO  Resistance  attributable  to  the  addition  of  gas  to  the 
mainstream  in  the  manifold  channel. 

R 1 1  Resistance  attributable  to  addition  of  gas  in  the  direc¬ 
tion  of  the  outlet  mainstream. 

R 12  Resistance  to  flow  through  the  outlet  manifold  chan¬ 
nel.  This  is  modelled  as  a  flow  through  a  pipe  with  a 
rough  surface. 

The  pressure  drop  across  each  resistance  is  related  to  the 
gas  flow  by 

\p=K\pv*  (1) 

where  K  is  the  resistance,  pthe  gas  density,  and  v  the  average 
gas  velocity.  In  some  cases,  the  friction  factor  can  be  deter¬ 
mined  from  formulae,  in  other  cases  it  is  necessary  to  rely  on 
tabulated  data  as  given  in  Ref.  (9].  These  data  are  likely  to 
be  measured  in  a  certain  set-up,  which  only  matches  the 
considered  configuration  to  a  limited  extent. 

The  flow  through  the  channels  adjacent  to  the  active  area 
will  be  laminar  in  most  cases.  Thus,  the  friction  factor  can  be 
given  in  analytical  form  by  Eq.  (2) 


K= 


LL 

Re  Dh 


(2) 


The  factor / depends  on  the  shape  of  the  cross  section  of  the 
channel  [  10],  for  a  round  channel /=64,  while  for  a  square 
channel  f- 56.8.  The  length  of  the  flow  path  is  /,  Re  is  the 
Reynolds  number,  and  Dh  is  the  hydraulic  diameter. 

These  formulae  were  used  to  evaluate  the  pressure  drop 
across  the  channels  (R6,  R7),  the  resistance  across  the  dis¬ 
tribution  area  (R4),  and  the  resistance  across  the  collection 
area  (R9)  in  the  plate. 

To  evaluate  the  resistance  from  the  divergence  (R2,  R3) 
and  the  convergence  (RIO,  R1 1 )  use  was  made  of  tabulated 
data  for  ducts  with  a  small  side  duct  that  extended  from  the 
main  duct  at  a  straight  angle.  Basically,  the  data  are  only  valid 
for  fully  developed  flow  in  the  feeding  main  duct,  which 
means  that  the  distance  between  subsequent  side  ducts  should 
be  at  least  5  to  10  main  duct  diameters.  In  the  case  of  a  fuel- 
cell  stack,  however,  this  distance  will  be  only  0.1  to  0.2 
diameters.  This  introduces  some  uncertainty,  which  can  be 
taken  away  by  a  comprehensive  analysis  of  this  configuration, 
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possibly  through  measurements  complemented  with  a  numer¬ 
ical  solution  of  the  energy,  impulse  and  mass  balances. 

The  same  holds  for  the  contraction  ( R5)  and  the  expansion 
(R8)  before  and  after  the  channels.  The  flow  is  likely  to  be 
fully  developed,  but  in  this  case  the  tabulated  data  are  valid 
for  one  diameter  change  only,  and  not  for  a  number  in  parallel. 

The  rough  pipe  flow  (Rl,  R12)  can  be  described  by  Eqs. 
(1)  and  (2)  for  low  Reynolds  numbers  (Re  <2300).  At 
higher  Reynolds  numbers  (Re >3500),  the  resistance 
becomes 

ff-t—  (3) 

where  £  is  the  friction  factor  presented  in,  for  instance.  Ref. 
Ill].  The  friction  factor  depends  on  the  roughness  of  the 
surface,  which  is  defined  as  the  ratio  between  the  average 
height  of  the  obstacles  in  the  flow  path  and  the  pipe  diameter. 
These  data  normally  do  not  incorporate  effects  from  splitting 
or  combining,  also  care  has  to  be  taken  when  interpreting  the 
results. 


3.  Numerical  approach 

The  model  solves  pressure  and  flow  distribution  for  the 
configurations  presented  in  Figs.  1  and  2.  Summarizing  the 
process  comes  down  to  solving  a  set  of  equations  of  the  form 

where  p  is  a  vector  containing  the  pressure  changes;  ^  is  a 
vector  containing  the  flows;  R#  is  a  matrix  containing  resis¬ 
tance  coefficients  multiplied  by  several  geometrical  and  phys¬ 
ical  parameters.  One  of  these  parameters  is  an  estimate  of  the 
flow,  which  is  identical  to  the  flow  in  the  elements  of  the  flow 
vector  <f>,  thus  the  problem  is  linearized. 

The  initial  flow  vector  is  found  by  assuming  that  all  plates 
get  the  same  flow.  This  produces  starting  values  for  the  coef- 


Fig.  2.  The  resistances  per  plate  are  added  up,  to  give  four  resistances  per 
plate.  They  ore  linked  together,  to  simulate  a  network  for  the  hydraulic 
resistances  in  a  fuel-cell  stack.  This  network  gives  rise  to  a  set  of  equations, 
which  after  solving,  gives  the  flow  distribution  in  the  stack. 


ficients  of  the  matrix,  RA,  and  for  the  coefficients  of  the 
pressure  vector,  p. 

Solving  the  set  of  equations  produces  a  new  flow  vector, 
and  subsequent  coefficients  for  the  matrix  and  the  pressure 
vector.  This  process  is  repeated  until  the  difference  in  length 
of  the  two  flow  vectors,  from  two  subsequent  iterations  rel¬ 
ative  to  the  length  of  the  flow  vector,  is  less  than  a  specified 
value.  The  matrix  contains  n  X  2n  elements,  where  n  is  the 
number  of  plates.  Because  it  is  a  band-matrix,  with  5  bands 
occupied,  only  5X2 n  elements  need  to  be  examined. 

The  flow  sheet  of  the  programme  is  described  in  Fig.  3.  To 
begin  with,  an  input  file  has  to  be  made  and  contain  data  that 
describes  the  geometry  of  the  separator  plate,  the  number  of 
plates,  the  physical  data  of  the  gases,  the  gas  flow  and  the 
flow  added  or  deducted  to  simulate  the  electrochemical 
reactions. 

To  initialize  the  calculations,  the  flow  is  distributed  evenly 
across  the  stack.  From  these  flows,  K  values  are  determined 
from  a  formula  or  a  table,  depending  on  the  application  for 
each  specific  case.  Then,  the  resistances  that  are  in-line  are 
added  up,  with  the  result  that  each  plate  now  consists  of  four 
resistances,  as  shown  in  Fig.  2.  As  part  of  the  linearization 
process,  each  of  these  resistances  is  multiplied  by  an  csti- 


Fig.  3.  Row  sheet  of  the  model. 


218 


RJ.  Boerxmu,  N.M.  Sammei  /  Journal  of  Power  Sources  63  (1996)  215-2J9 


Tabic  1 

Viscosit)',  density,  composition  and  flow  of  the  gases  at  in-  and  outlet  of  the  stack.  The  properties  at  the  cathode  outlet  differ  for  each  load  setting,  but  the 
difference  with  the  inlet  properties  hardly  affects  the  results 


Inlet  anode 

Outlet  anode 

Inlet  cathode 

Outlet  cathode 

Viscosity  { Pa  s) 

21.TXI0"4 

48.1  Xl0_h 

46.8  X  10”  A 

depends  on  flow  and  current 

Density  (kg/nv  ) 

0.022 

0.21 

0.31 

depends  on  flow  and  current 

Composition 

pure  H2 

20%  H2, 80%  H20 

air 

depends  on  flow  and  current 

Row 

set  for  80%  utilization 

set  for  80%  utilization 

set  for  A  T=  150  °C 

set  for  AT  =  150  °C 

mated  value  for  the  flow.  This  produces  the  5x2 n  band- 
matrix  mentioned  above.  With  the  estimate  for  the  flow,  the 
pressure  drop  for  each  element  can  be  calculated,  which  gives 
the  left-hand  side  of  Eq.  (4). 

The  set  of  equations  is  solved  with  a  standard  routine  for 
band-matrices,  which  gives  the  new  value  for  the  flow  vector. 
The  process  of  filling  the  resistance  matrix,  and  the  pressure 
vector,  and  the  subsequent  solving  of  the  flow  vector  is 
repeated  until  the  convergence  criteria  is  reached. 

Convergence  was  significantly  improved  by  applying 
some  under-relaxation  to  the  flow  vector;  the  average  of  the 
new  and  old  flow  vectors  was  taken  as  the  new  flow  vector. 

The  programme  was  written  in  Turbo  Pascal;  the  program 
code  requires  107  kB  of  memory  for  execution  and  10  kB  for 
data.  Large  arrays  are  allocated  dynamically  through  the  use 
of  pointers  and  allows  data  to  be  stored  for  300  plates.  The 
data  occupy  approximately  640  kB.  Convergence  is  usually 
reached  within  seconds. 


4.  Case  study 

The  analytical  model  that  was  discussed  in  Ref.  |8]  has 
been  used  to  study  the  distribution  of  gas  flow  in  the  SOFC 
stack  of  Tokyo  Gas  [  12,13] .  Details  of  the  design  were  not 
available  for  this  study,  therefore  an  assumed  configuration 
was  used.  Because  the  numerical  model  allows  for  more  detail 
than  the  analytical  model  described  in  Ref.  [8],  additional 
assumptions  had  to  be  made  regarding  the  geometry  of  the 
stack,  the  in-  and  outlet  gas  compositions,  the  viscosity  and 
the  density.  Table  1  shows  the  gases  and  their  properties  used 
for  the  calculations.  7  he  inlet  data  are  for  850  °C,  the  outlet 
data  are  for  1 000  °C. Tie  cathode  flow  is  assumed  to  remove 
all  the  process  heat,  and  its  flow  is  set  accordingly.  The  fuel 
utilization  is  assumed  to  be  80%  and  is  kept  constant  by 
adjusting  :hc  flow. 

At  maxir.um  power  output,  the  cell  voltage  is  0.5  V,  giving 
1  A  cm-2  1 13],  which  corresponds  to  0.5  W  cm-2.  At  this 
setting,  the  heat  production  is  at  a  maximum  at  0.78  W  cm  "  2. 
For  this  setting,  the  cathode  and  a  node  gas  distribution  were 
calculated  and  arc  presented  in  Figs.  4  and  5.  together  with 
results  for  80  and  55%  of  the  maximum  power  output. 

It  can  be  concluded  that  the  flow  distribution  for  the  cath¬ 
ode  is  most  likely  to  be  the  problem,  specifically  at  maximum 
power  output.  When  the  output  was  reduced  to  80%  of  the 
maximum,  however,  the  flow  distribution  was  greatly 


improved.  The  flow  distribution  at  the  anode  side  is  not  likely 
to  be  a  problem. 

With  the  analytical  model  [8],  it  was  calculated  that  for 
the  maximum  output,  the  ratio  between  the  flow  through  the 
top  cell  and  the  average  cell  would  be  0.89  for  the  cathode. 
The  computer  model  shows  that  this  ratio  is  about  0.85.  This 
implies  that  the  top  cell  gets  a  15%  lower  flow  than  the 
average,  and  thus  z  15%  higher  temperature  difference. 

The  difference  is  caused  mainly  by  the  fact  that  the  ana¬ 
lytical  model  uses  an  average  viscosity  and  density,  whereas 


-100%,  Intel 
-100%.  outlot 
-80%,  Intel 
—80%.  outlet 

—  55%,  Intel 

—  65%.oullol 


cell  number 

Fig.  4.  Gas-flow  distribution  at  the  anode  side  under  various  loads.  The  flow 
is  set  for  80%  utilization.  The  flow  increase  at  the  outlet  compared  with  the 
inlet  results  from  the  density  decrease,  which  results  from  the  tempetature 
increase,  due  to  the  electrochemical  reactions. 


-100%,  Inlet 
-100%,  outlet 

-80%.  Intel 
-W~80%,  outlet 
-6P%.  Inlet 
-55%.qullct 


cell  number 

Fig.  S.  Gas-flow  distribution  at  the  cathode  side  under  various  loads.  The 
flow  ts  set  for  an  average  temperature  difference  of  1 50  °C.  The  flow  increase 
at  the  outlet  compared  with  the  inlet  results  from  the  density  decrease, 
resulting  from  the  temperature  increase  and  the  molar  flow  increase,  resulting 
from  the  electrochemistry.  The  net  result  is  an  increase  of  the  volume  flow. 
The  high  flow  at  the  lower  end  of  the  stack  will  effect  a  lower  temperature 
difference  between  in-  and  outlet  than  at  the  high  end.  This  will  cause  heat 
transfer  from  the  top  of  the  stack  to  the  bottom. 
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the  computer  model  uses  different  values  for  in-  and  outlet. 
When  average  values  are  used,  the  computer  mode!  and  the 
analytical  model  produce  practically  the  same  results. 


5.  Future  work 

Clearly,  more  reliable  results  can  be  obtained  by  modelling 
the  flow  and  heat-transfer  phenomena  in  a  dedicated  com¬ 
putational  fluid  dynamics  (CFD)  programme,  such  as  Fluent 
or  Phoenics.  These  programmes  arc  capable  of  simultane¬ 
ously  solving  mass,  impulse  and  energy  balances  for  practi¬ 
cally  any  configuration.  They  can  certainly  provide  more 
accurate  data  on  items  such  as  splitting  or  combining  a  flow, 
which  would  add  to  the  reliability  of  the  model  presented  in 
this  paper.  Alternatively,  if  one  is  successful  in  combining  a 
fuel-cell  model,  with  a  CFD  code,  then  a  very  accurate  stack 
model  may  be  obtained.  For  instance,  heat  transfer  along  the 
height  of  a  stack  can  be  modelled,  something  which  has  not 
had  much  attention  until  now,  but  is  justified  because  of  the 
unevenness  of  the  flow  distribution  that  can  be  expected. 

Also,  experiments  may  provide  data  that  can  he  correlated 
to  numerical  findings.  The  set-up  for  these  experiments  is  not 
necessarily  complicated,  as  they  can  be  carried  out  at  room 
temperature,  using  air  as  the  flow  medium,  and  built  from 
plastic  components.  The  inferred  friction  factors  should  also 
apply  at  high  temperature  and  using  hydrogen  a?  a  fuel. 


6.  Conclusions 

I .  A  computational  model  was  developed,  describing  the 
gas  flow  and  pressure  relationship  in  terms  of  resistances, 
which  is  numerically  solved  by  an  iterative  process. 


2.  In  the  design  of  internally-manifolded  SOFC  stacks,  the 
distribution  of  the  cathode  gas  needs  more  attention  than  the 
distribution  of  the  anode  gas. 

3.  Results  from  the  numerical  and  analytical  model  match 
well.  This  gives  confidence  in  the  reliability  of  both  models. 
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